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STABILITY ANALYSIS FOR LAMINAR FLOW CONTROL - PART I 


David J. Benney and Steven A. Orszag 
Cambridge Hydrodynamics / Inc. 
Cambridge, Massachusetts 02139 


1. INTRODUCTION 


The theory of stability of quasi-parallel flows has 
been subject to extensive study during recent years. The 
problem is nonlinear and analytical progress is difficult 
without the introduction of some form of approximation. 

For this reason, much attention has been given to the linearized 
problem so that the study of Or r- Sommer fe Id- like eauations 
has been the dominant theme. In particular, accurate methods 
for the calculation of eigenvalues and eigenfunctions has been 
and remains an important task. One of the two main goals 
of the present work is to assess and recommend efficient 
and accurate numerical methods for the calculation of 
eigenvalues and, when necessary, eigenfunctions. A full 
discussion of these matters is given in Sects. 5-8, 

Beyond the linear regime, various nonlinear theories 

have been proposed and each of these have advantages and 

defects depending on the particular physical problem under 

study. Here we are concerned with boundary layer stability 

predictions and for this purpose the technique which will 

be adopted in Sects. 2-4 is the standard one in which 

the nonlinear evolution of the basic amplitude is studied 

by simple perturbation techniques. The idea of this method 

1 2 

originated in the work of Stuart and Watson • When 

the non-parallel aspects of the flow are incorporated into 

this system of equations, the result is a nonlinecir space— 

3—8 

time evolution equation for the amplitude. 


In Sects. 2-4, we formulate a theory for the nonlinear, 
nonparallel stability of a very general class of boundary 
layer flows. The effects of nonlinearity, compressibility, 
and the quasi-parallel nature of the flow are all included. 
The theory given here is an important first step in trying to 
obtain a realistic approach to the prediction of transition 
in three-dimensional boundary-layer flows like those 
encountered on laminar- flow- control aircraft. 


2. STABILITY OF QUASI-PARALIiEL INCOMPRESSIBLE FLOW 


It is instructive to work out the stability theory 
for incompressible flows before proceeding iii Sects. 3-4 
to the more complicated case of compressible flows. Here 
we formulate the nonlinear, nonparallel stability theory 
of two-dimensional disturbances in incompressible boundary 
layer flows. The basic ideas used here carry over directly 
to the problem of three-dimensional disturbances in a 
compressible flow. 


It is convenient to use a s treamf unction i}/ , so that 

the equation for i|; is 

If 5(y/X) is the basic flow with X as the slow variable 
(to indicate non-parallel effects), and t|^(y,9fX,T) is the 
perturbation s treamf unction (with T as the slow time variable 
and 0 as the phase of the perturbation) , then 

(0 iS ■ + 0 ) (iIj ^ + 0^ - 0 ip 

^ x^y T ^yyB x ^006' x^yyy ^0 

2 4 

^^yyyy x ^yy00 x ^ 0000 ' 

+ u[20„0 ^ ™ + 20 0 

^ T x^00x yyT X ^00T T XX X XT 69 

+ 5 {tl; + 30^ + 30 0 +4^ , 

y ^yyx x 09x x xx^00 yyx y 

~ ^x^^yyy ^ ^x ^y0 9^ "" ^^^®x^yy0x ^ ^^xx^yyB^^ 

* 'IVyWyyeT^ l'eee> ' ®x»6<*yyy * ®x »yee” ■ " 


( 2 . 1 ) 


( 2 . 2 ) 
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II 


Here we have assTimed that the phase 0 is of the form 

0 * 0{X,T)/u , 

X = ]ix, \x = 6/L where 6 is the boundary layer thickness and 
L is a measure of the distance over which non-parallel flow 
effects are important, and e is a parameter specifying the 
amplitude of the nonlinearity. In formulating (2,2), we have 
also used the relations 


A = + 0 -L A 

^3X 90 ' 9t 


= 9. 


=*T 39 


In order to solve (2.2), jp is expanded as 

+ etAA%(°) + A2^<2)^2i9 ^ ^*2^(2) V2i9j 
+ e2lA3^(3)^3i0 ^ ^*3^(3) V3ie ^ ^2^*^(U)^i0 


.(i).ie 


(D* ^-ie ^ 


+ 'e^'' + A^ e + Aij; 


+ A ii«2 


e-i®] + . . . 


(2.3) 


where A(X,T) satisfies the slow space-time scale equation 

At = Oj^Ajj + a^A + XA^A* 


(2.4) 


If nonlinear and nonparallel effects appear at the same order, 

2 

we must take y = e . Making this assumption, the rest of 
the calculations proceed by equating terms order by order 
in (2.2) . 


4 



The first few boundary-value problems obtained from 


(2.2) are: 


i(0 il + 0„) - 0^ - i0 5 

' x^y T' ^yy x ^ ^ x^yyy^ 


- V(*I« - ♦«> + e‘ ♦< 1 >) . 0 

yyyy x ^yy x ^ ' 


(2.5) 


^yyyy ^yy x ^ ^ '^^yy 


®x »‘^’> 


+ i»(l) *(»<!) - g2 ^(1>, 

yyy x ^y ' 


- - 0%(1)*)] = 
^yyy x ^y ^ 


( 2 . 6 ) 


- - 80^ + 160^ 

^ xrrr^rrr v ^ ^ ~ J 


yyyy 


X 'yy 


+ 10 [Ip (1/; - 0^ - 0^ 

X y yy X ' yyy x y 


= 0 


(2.7) 


H- 0^) (^yy 


- 02 ^(11)) - i0 5; ^(11) 

X ^ ' x^yyy^ 


- ^<Cy - 

+ e^<>(-i)<>* - 8^ ♦"'*) * ♦f>*(2i)(*^> - 4el 


- - e? *W*) .i» 


ryyy 


X ’^y 


-'•(^>*(ip(2) . 40^ ^(2) 


ryyy 


X ^y 


^ ♦<<”(i)wiV - e! ♦<^>) - i*<'>> ♦'W) - 


ryy 


yyy 


( 2 . 8 ) 


(2.9) 




X'^yyy'^l 

lyyyy x ^lyy x ' i yy x ^ ' 

- + ?y{<l'yy^- 3 O ^ ^ ^ ^ ~ V4i0^</^<j) = 0 


i(0 ^ + 0„) (<()i^^ - 0^ Oii^^) - ie iS 

' x^y T' '^2yy x ^2 ' x^yy^2 


- - 20J + ©4 ^<U) + a-(,i»iV - 0 ? 


2yyyy 


( 1 ) 


X '^2yy ^ ^x '*'2 

( 1 ) 


2''*'yy 


- 20„0 i|)' '(0„0 + 20 0 )i(i 

T x^x T XX X TT 


+ - 30^ - 30 0 

y yyx x ^x x xx^ ^yyx^y 

“ U; - 0^ - iv(40 + 20 

^x'^yyy x, '^'y ^ x^yyx ^ ^ xx^yy ^ 

+ x0 $ • 0^ - i0 f = 0 

x^y'^yy x ^ ' x^yyy^ 


Eq. (2.5) is the homogeneous problem 
^ 0 


where jC is the Orr-Sommerfeld fourth-order differential 
operator. The usual form of the Orr-Soiranerfeld equation 
follows by making the identifications 

k = 0^ , 0) = - 0^ . 

The boundary conditions on are that =0 

at y = 0,» . When (2.11) is solved subject to these 


( 2 . 10 ) 


( 2 . 11 ) 



boundary conditions, there results an eigenvalue condition 
of the form 

f(0^,0^,X) = 0 (2.12) 

that must be satisfied for a nontrivial solution to (2,11) 
to exist. At each position X , we can calculate 0^ given 
0^ or vice versa from (2.12) . If 0^ is assumed given and 
real, then (2.12) yields a temporal stability analysis; if 0^ 
is assumed given and real, then a spatial stability analysis 
results. In either case, it is assumed that conditions are 
close to neutral so that any weak amplification or damping 
can be incorporated later into the evolution equations. 

Eqs. (2.6-7) are inhomogeneous problems that can be 
solved at each location X. The more interesting problem 
arises with Eqs. (2.8-10) since they are each of the form 

= K F + G (2.13) 

Here K(X) is either X, a^, or U 2 in (2.4), and F and 
G are known functions. In order to determine K, we invoke 
the Fredholm alternative so that 


K 


/G X dy 
/FX dy 


(2.14) 


where X is the solution of the homogeneous adjoint problem 

7)lX =0 (2.15) 

That Ls, Tfl is the adjoint of ^ and the boundary conditions 

for % are adjoint to those for . 

By this technique, all the functions appearing in (2.4) 
are found and hence the evolution of the amplitude A can be 


detesrmined. 



3. yHREE-DIMENSIONAL NONLIIOEAR NONPARALIiEI. STABILITY THEORY 
OF COMPKESSIBLE BOUNDARY EAYER FLOWS 

The equations governing the motion of a compressible 
fluid are as follows: 


P 


DE 

Dt 


It * slj ‘‘’“i’ 


Du. STj. 

-dT - *-55^ 


» - If J t 5I- 


p = F(p,T) 


(3.1) 

(3.2) 

(3.3) 

(3.4) 


These equations are those of continuity, momentum, energy, 
and state, respectively. The notation is standard, viz. 
p, P/ T, V = (u2^,U2/U2) / X = {X 2 /X 2 /X 2 ) , t, E, k, X = (X2#X2,X2) / 
T^j, and $ denote pressure, density, velocity, spatial 
coordinates, time, internal energy, coefficient of thermal 
conductivity, external force, stress tensor, and dissipation. 

In addition, we use the following subsidiary relationships: 


5 3u, 9u . 9u . 

T.j = - (p + ^v ^)6.. + v(^+.^) 


(3.5) 


E = C^T 


$ 


e. .e. . 
2 X3 ID 



(3.6) 


(3.7) 



where 


9u. 3u . 
®ij ” 3x7 33 c^ 


(3.8) 


A == div V 


(3.9) 


Here v and are the coefficient of viscosity and the 
specific heat at constant volume. It is further assumed 
that the three transfer coefficients k, v, and c^ may be 
assumed functions of temperature T alone. 

Armed with these basic equations, the aim is to develop 
a theory for the nonlinear stability of a basic three-dimensional 
boundary- layer flow. The boundary will be taken to be at 
y = 0 so that the basic flow will have slow variations in 
the X = and z = directions, together with possible 
slow temporal variations. This feature is most easily 
accounted for by the introduction of slow variables X, Z, and 
T* (do not confuse with T) where 

(X,Z,T*) = y(x,z,t) / y << 1 

The small parameter y is the usual one on which boundary layer 
theory is based. 

The unperturbed boundary layer solution will be denoted 
by overbars. This zeroth order solution is taken as given. 

In addition, the first-order corrections to the basic state 
do influence the evolution of disturbances - 

The basic state is then given by the velocity components 


9 



V = (u(y,X,Z,T') , yv(y,X,Z,T’) , w(y,X,Z,T')) (3.10) 

and density, temperature, and pressure fields 


p(y,X,Z,T')» T(y,X,Z,T'), p(y,X,Z,T') 


(3.11) 


For the stability analysis, each function f is replaced 
by f + f and perturbation equations are then derived- Third- 
order terms in amplitude must be retained in the present theory. 
To be specific, we perform some preliminary calculations. 

For example, since 


then 


and therefore 


E = c^(T)T, 
E = c^(T)T, 


(3.12) 

(3.13) 


E = c^(T + T) (T + T) - C^(T)T . (3.14) 

Taylor expansion of this latter result gives 

E = + E^T^ + ... (3.15) 


where 


E 3 = cy2 +c^"T/ 6 . 

Here primes denote derivatives of c^(T) = c^. 

In the equation of state for the basic flow 


(3.16) 


p = F(p,T) , 


(3.17) 


10 



so that the perturbation pressure is given by 


p = F(p+p, T+T) - F{p,T). 

Taylor expanding, it follows that 

P = + |[F,nP^ + FnpT + 


01 " " 2 ^" 20 ^ 




Where F „ = 3 "‘‘‘'"f ( p,T) / 3 p™ t” 
mn 

In a similar way, the stress tensor is expanded as 

9 _ - au. _ „2 3 u. 

■'ll = -■'ijlP ^ l‘<“ 3 ;;^* li5^> * "‘■’^3;^* "" 3-33q;> 

+ (U». 10 ^+ u"’ J 3 

+ 2 3Xj^ 6 3Xj^' ^ ^ 


8 u. 3u . 9u. 9u . 




3u. 3u. -„ - 3u. 3u. 

* <>■'’'<357 * 3;^> * ¥ <33q- * 35 tJ>> 


_ 2 3u. 3u. 3 Su^^ 3u. 

+ ^ < 3 ^ 33^> + V ■' <3^ + 33jJ> > ' 


and the dissipation function i is expanded as 


(3.18) 


(3.19) 


(3.20) 



* “ v) (div v) ] + u'T[| i^g - |(div v)^]} 


2 a3 3 ’ 


+ e^g - j(div V)^l + V>'T[ej^ge^g - ^(div v) (div v) ] 


+ ^ l| -«3 


+ {^'T[| - |(div v)2] + [i„ge^g - | div v div v] 


+ ^ t|«a3 . 


(3.21) 


Other functions which occur in the fundamental equations 
can be expanded in this way and need not be recorded explicitly 
at this stage. Of more importance to our subsequent analysis 
is the fact that perturbation quantities, say g, will have 
both fast and slow space-time scale variations (see Sec. 2) . 

In order to incorporate this feature, it is again desirable 
to introduce a phase function 


0 


0(X,Z,T') 


(3.22) 


so the perturbed quantity g is expressed in the form 


g(y,e,x,z,T') 


(3.23) 



It follows that derivatives of g depend on both the slow 
and fast variables* For example. 


If = ©xge + <3*24) 

M = ®X 500 + (20^93^ + Gjjjjgg) + (3.25) 

oX 

Other derivatives can be calculated in a similar manner. 

Each perturbation function is expanded and the nonlinear 

parameter e is introduced to order products of perturbations. 

The two small pareimeters p and e are used to order the 

non-parallel and nonlinear effects, respectively. In order 

2 

for these effects to be in balance, the choice y = e turns 
out to be the appropriate one. 

The equations are now rewritten in this notation and 
only the terms consistent to the order of the calculation 
are retained. 

For the continuity equation, we obtain 


13 



, + wlp^ + (pu + pu)jj + (pv)y + (pw + pw)^3 

+ e[0j^(pu)Q + (pv)y + 0j^(pw)g] + ep[(pu)^ + (pw)^] = 0 (3.26) 


For the X momentum equation. 


PtQ^, + (0jjU + 0zW)]Ug + p ly V 


r" *9 *“9 “”9 *9. ” 9u ~ 9u 


+ (u^ + uUjj + vUy + vni^)p] 


^ du. ** ^ 

+ e[p(0jjUUg + Wy + O^'^q) + ®tP'^0 ■*■ 3y '*' ®z'^‘^'^0^ 


+ e'^TpCOj^u + 02 W)Uq + pvu ] 


l-e^Pg - yPx “ I ^ 0 x^®x '^0 0 + ®z''ee + '^ye^ 


I ^^f^®x" 0 x ■*■ ®z'' 0 x ■" ^xx'^e ®xz"e^ 


I v*pT^(0^Ug + 0 ^w^ + Vy) + v202 Uqq 


+ yv(2) (0^ug^.+ 0^^ug) + y2VT^0^UQ 


14 



+ Pi -I V(Ug^ + Wg^)e^ - I V'Tg(5^ + ) 9^ 

+ 2ve^"ex + 2v*e^5^Tg] 

+ el-| V*0j^{T(e^Ug + e^Wg + Vy))g + 2V*0^ (TUg)g] 

+ e^I-| v"9x^^^(®x"e ®z''e * v"0^(T^Ug)g] 

- 

+ lV(Uyy + 0^Vyg) + V'Ty(Uy + 9 ^ V g ) V ^ T 

+ V'UyTy + V"TyUyT] 




+ e {v'T(Uy + 0j^Vg) + \ 




+ V 9jj(0 jjWgg+ ©2 Ugg ) 


+ y [V (©^Wg^ + O^Ug^ + -02 3 , Wg +©22 Uq) + V T ^ ^ Wg+- 0„U g ) 


+ V (W^g + U^g) ©2 + V 0^ Tg(U^+Wj^)] 


+ e© 2 lv’ TC©^ Vg+ £> 2 Ug)]g 


+ e^©j,[v"-|- (0'x''e'^®z 


(3.27) 
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The y momentum perturbation equation is 


P[©J + + ©2^)]Vg 


+ U[p(v^ + ^ + V ^ + w j^)v + P V] 


+ etp((0^u + 02W)Vg + Wy) + 0 ,j.PVq] 


+ e^[p((u0^ + w9^)Vg + w )] 


- Py - I ^(®x"ey + ®z''0y + ^yy^ 


- I ^’Ty^Ve + ®z''e + ^y^ + 2 v'TyVy] 


+ p[ v(u„ + W„) - -I- V' (u^ + v„ + W„)T + 2V'v„T] 


3 '''“x " "z' 3 '“x ■ ’y ■ ”z' 


y y 


+ e[ - j v'T(0^Ug + 0j^Wg + Vy) + 2v'TVy]y 


+ e^I-j v"T^(0j^Ug + 0jjWg + Vy) + v"T^Vy]y 


+ •" ®z^e 6 ^ ''’Ve’ 


+ V[V(Wy^ + 0^Vg^ + 0^^Vg) + VT^(Wy + 0^Vg) 


+ ^'WyTj, + v'WyjjT + v«T^WyT + ©2 ^Vq ^3 


16 



+ e 0 ^[v-T(Wy + 0^Vg) + ^ WyT^le 


+ Qz'^q) + ^ V ^]0 


+ + e^Vgg) + V'GyTg] 


+ p[v(Uy^ + 0^Vg^ + 0 ^^Vq) + v'T^(Uy + 0^Vg) 


+ VUyT^ + v'Sy^T + v'T^SyT + 0 ^vVq^] 


+ e0^[vT(Uy + 0^vg) + ^UyT^lg 


2 — — 

+ £ 0 ^[v" (U„ + 0 V.) + - 2 - U T 


■'x'-' 2 '“y ''x’e 


y J e 


u . 


( 3 . 28 ) 
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The 2 momentiim equation is analogous to the 
corresponding x equation and takes the form 


PiQrp + + P Wy V 


+ V[p( ^+GA + v^+w ^)w + p|g„+p||u 


+ p(w^ + UW^ + VWy + ww^)] 


+ eip(0jjUWg + vwy + ©^w^e) + ®TP''e ■*■ a3^ 


+ e [p(0jjU + Q^w)Wq + pvWy] 




+ ^t^^Qz'^ex + ^x'^ex ^zx’^e ^xx'^e’ ■*• ■" ®x"e^ 


+ ^^"z 0 + ''xe^^x + ^’®x'^e^^z + ''x^^ 


+ e0^[VT(0^U0 + 0xW0)]e 


+ e\l ^ ( 0 ^Ug + 0 ^Wg)]g 


+ V(Wyy + 0^Vgy) + VTy(Wy + Q ^ 


+ v"TyWyT + v'WyyT 
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+ e[V'T(Wy + e^Vg) + — V^ly 

*> *> V”' — 3, 

+ e^l-^ T2(Wy + e^Vg) + -g- WyT ]y 

+ I-Pe®2 - I ^®z(Vee ^ ®z''ee Wgg] 

+ Pt-P^ - I V(0^Ug^ + 0^Wg^ + Vy^ + 0^^Ug + 0^^Wg) 

" I ^’^2<®x"0 •" ®z''0 ■" ■" ^^©zWg^ + (W^)y 

+ 2V0^^Wg + 2V-T^0^Wg 

- I ®Z^("0X ■" ''0Z^ ' I ^*®z'^0("x ^y "■ 

+ ®z2^"ez 2v'w^0^Tg] 

+ V'T(0^Ug + 0^Wg + Vy) + 2V'0^TWg]g 

+ c\l-J V"T2(0^Ug + 0^Wg + Vy) + V"0^T2wg]g . 


(3.29) 


19 



The energy equation is 


^®z^0 


+ V[pCy{T^ + uTj^ + vTy + wT^ + nT^ + wT^}] 


+ etpc^Ce^uTg + vTy + © 2 '^ 0 > 


+ (c^P + pc^T) (0 ^Tq + u 0 j^Tq + wG^Tg + TyV) ] 


+ e^[(c^P + PC^T) (©x“^e ■*■ ■*■ ®z^^ 0 ^ 


r -^v „2 


+ (C^PT + p -f- r") (G^Tg + Gj^uTg + G^wTg + TyV)] 


v[2Sy(Uy + Gj^Vg) + 2Wy(v7y + B^Vg)] + ^'(11^ + w^)! 


+ p[VI25yV^ + 2 WyV^ + 2 (u^ + w^) (G^Ug + G^Wg) 


+ j { 2u^0^Ug + 2VyVy + 2w^G^Wg 


- CS^Vy +.VyG^Ugl ^ (w^Vy + VyG^Wg)}J 


+ e[V{(Uy + OjjVg)^ + (Wy + e^Vg)^ + (Gj^Ug + G^Wg)2} 


* T '4 + ''J * “e - ‘®x“e * w^'y - V2"e”8> 


+ v'T{2Gy(Uy + 9jjVj) + 2Wy(Wy + OjVj)> + iy T^(Uy + Wy)l 


20 



+ e^[v’T{{Uy + + (Wy + © 2 ^ 0 )^ + 


+ i V'TCej u2 + vj + ©2 „2 _ (e^Ug + e^Wg)Vy - Oj^e^Ug^g}] 


- 2 

+ {2Uy(Uy + e^Vg) + 2Wy(Wy + 02V0)} 


V"' r-2 . -2, 

+ - 6 - 


- TFg^(0^Ug + Vy + 0^Wg) 


+ y[ - T(u^ + Vy + W^) + Fq2T) - T(u^ + w^)Fgj^ 


- ’’^01 ''z>J 


+ £[ - T(0^Ug + 0^Wg + Vy) (PF^^ + TFg2> " T(0^Ug ^ 0^Wg + Vy ) F g ] 

• 2 2 

+ e2[_T(0^Ug + 0^Wg + Vy){ ^ F2^ + PTF 22 + Fq 3 5^ ) 


- T(0^Ug + 0^Wg + Vy) (PF^J^ + TFq2)] 




+ PKkG^Tg)^ + 0 ^kT^g + k-T^ 0 ^Tg + ( kO ^ T g ) ^ 




+ e[k'ej <TTq)q + ic-e^ (TTg)g + (k'TTy)^] 


+ e2[^(02 + e2){T2Te)g+ (H:^)yyl . 


Finally, the equation of state is 


= ^’lO + ^Ol'" + ^ ^ 20 P^ + + I 

+ ^ ^3QP^ + 1 I ^12^’’^ + I ^03^^’ • 

The perturbation equations (3.26-31) are to be solved 
subject to appropriate boundary conditions. At y = « , it 
is clear that all perturbations should decay to 0 so we 
require 

0 

(p,u,v,w,p,T) '►0 (y -► «) 

At the rigid boundary y=0, u=v=w=0 and the 

temperature condition will be of the form 


oT + 




0 


The case B = 0 corresponds to an isothermal wall and the 
case a = 0 corresponds to an adiabatic wall. 


(3.30) 


(3.31) 



I 


4. FINAi; EQUATIONS FOR THREE-DIMENSIONAL NONLINEAR NONPARRLT.T!!T. 
STABII.ITY THEORY OF COMPRESSIBLE BOUNDARY LAYER FLOWS 

First consider the linear problem so that e *= u =» 0 . 

We ask for solutions in which each amplitude function is 
written in the form 

h(0,y,X,Z,T‘) = Ah^^^ (y,X,Z,T*)6^® (4.1) 

where the envelope conplitude is to be determined as a function 
of the slow variables, i.e., 

A = A(X,Z,T') (4.2) 


The respective equations are: 

1(0^, + U0JJ + W02)p^^^ + i0jj + 102 pw^^^ 




(4.3) 


ip(0T* ®x“ 




+ Y0„(0„u^^^ + 0^w^-^') - {y(u'‘/' + i0yV'-^'} 


'Z- --y 

,(1). _ 


( 1 ) 


z' z 


-Y 


'X 


- (Y' UyT^^^) 


= 0 


(4.4) 
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0 


+ 02W)v^^^ + 

- i(0jjV'Uy + 0j,v'Wy)T^^> - + iOgV^^’) 

- + i02W^^^)J = 

ip(0^i + 0^u + 02W)w^^^ + 

- [j(w^^’ + i02V<^^)jy + 2 - (V'WyT<^^^ 

- + 0jjW^^^ = 

PCy|i(©^i + U0JJ + w©2)T^^^ + 

- y{2Uy(u^^> + + 2wy(w^^> + i02V^^>)} 

“ Y*(u| + wj) - Fqj^(10jjU^^> + + x0j,w^^^T 

- (IcT^^^y + l?(0x + ©2)T^^^ = 0 


(4.5) 


0 

(4.6) 


( 1 ) 

(4.7) 
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I 


p<» - P,„P«> - - 0 ,. 


The above set of equations (4.3-8) constitutes the standard 
eigenvalue problem for determination of the linear stability 
of a given boundary layer configuration. Rather than eliminate 
variables it is convenient to consider this system of equations 
to be written in matrix form 


(A, + C 

-1 -1 dy -I'ry. 


where ) and the 6x6 

matrices are readily written down from the pre- 

vious equation. 

The solution of the linear eigenvalue problem gives a 
relationship between -0_, r each position 

and time X,Z,T' . Of course the quantities ” 

0J, are the local frequency and X and Z wave numbers. At 
this stage we make no distinction between temporal and spatial 
amplification. In actual computation where the initial insta- 
bility is being followed and we are close to neutral conditions 
the variable 0 is treated to be approximately real. Note 
that at this stage of the calculation the wave packet amplitude 
A(X,Z,T*) is artJitrary. 

In order to proceed into the nonlinear problem it is 
necessary to return to (4.1) and replace this equation for any 
perturbation function by an expansion of the form 


1 



= Ah«>e“ t A*h‘»V“ 


t ♦AA-h"*') 

+ u[j;Ajjh|^’ -t Ajh‘^’ ■!■ Ah^*’)e''® 

* ,A5h<»*.A*h<»*.A.h<^>*,e-^«] 

+ c2[A2A.h'l»e« + A*a\<“>%-“ 

+ A\'=‘>e“« + A.\”>%-^“] 


2 3 

+ 0(epry ,e ) 

This expansion for any function h can be replaced by 
the identical expansion for the vector 

H = (p,.u,v,w,T,p) , 

'Vr 

where is defined in (4,9). 

The expansion and truncation procedure is straight- 

2 

forward and based on the fact that p and e are of the 
same order. The process necessitates that a corresponding 
expansion be invoked for A(X,Z,T’) of the form 

yA^, = y(a^A^ + «2^z ^ e^XA^A* , 

where r U 2 r »3 r and X are scalar functions of 

X, Z, and T* . These functions in this evolution equation 


(4.10) 


(4.11) 


(4.12) 
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are determined by orthogonality conditions associated with 
the other boundary value problems. In particular those for 

' are .connected to these 

functions and in addition the non-resonant boundary value 
problems for and are needed. 

( 2) 

First we write down the inhomogeneous problem for H , 
namely 

Here the solution used is that 


- 0 


( 4 , 13 ) 


IS 


( 2 ) _ 


(p(2), v<2), w<2). 


T<2), p(2)) 


(4.14) 


while the matrices A^, Bor and Co are identical to those 

f\^ r\j-^ f\j Z 

for A^r B^r and except that the variable 0 is every- 

(2) 

where replaced by 20 . The vector ^ is readily found 

and if 


k(2) 

At 


( 2 ) .( 2 ) .( 2 ) .( 2 ) .( 2 ) ( 2 ). 

, JC3 / , Kg , JCg ) , 


(4.15) 


we find that 


:f> = 2ie^p«>u<W . 2iv‘"’w‘'’ * (P<»v'»), , ,4.16) 
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.< 2 ) 




y * 


+ a (e^ + U0JJ + W0g)u<^^ + UyV<^^) 


] 


+ jY' 0x + i02W^^^) 


2y’0^ 


JT ^ J; jf 


2iY'O2T^^^i(02U^^^ + ) . 


(4.17) 




(1)„(1)1 


2 

- ‘2i0jj{Y'T^^^ (Uy^^ + i0j.V^^^) +jY"UyT^^^ } 


+ {|Y*’r<">(i0^u<l) ^ _ 2Y.T<">v^">}y 


- 2i02{Y'T^^^ + i02V^^^) + |Y”WyT^^^ } (4.18) 


28 



(2) 


= |p(i0jjU^^^w 


(1) + ^ ie „(i)^ 

y ^ 


+ i(0^ + U0JJ + w02)w^^^ + 


2i0jjY'T^^^i(02U^^^ +.0jjW^^^). 


+ ie + I y"w } 

y - " ^ j j 


2i02 <- I-y’T^^^ (i0jjU^^^ + ) 


+ 2Y*e2T^^^iw^^^ }_ , 


(4.19) 


<2) = [^y(i0^u<l>T(l> + * i0,w(i)T.(i>) 


+ (C„p^^^ + pc T^^^)(i(0_ + U0„ + w0_)T^^^ + T 


- Yt(u<"> + i0^v(i>)2 ^ ^ iO 

y X y ^ 


- (6ju'i> . e^w'*’)' 


t4(-e^'“^ .„<i>^..e^„a)^ 


- i(e,u'^> . e,w'^')vy' . e,e,u>*'w'*’ ) 

A ^ y 'A It. 


- 2y'T^^^ {u„(u^,^^ + + w„(w^^^ + i0„v^^^} 


y y 


y y 


- i . wj) 




(4.20) 




- K'i(eJ + ©2)T^^^ - k' , 


- ”2*^20*^ ^ * 2*02^ 


(4 


With these specifications and the appropriate boundary 
conditions the inhomogeneous problem defined by (4.13) can be 
solved. 

In the Same way 


. <3“ <3 + c - 

^^0 ^ ^0 dy ^ ^ 


= 0 


(4 


is the appropritate 

equation governing 

mean motion* Here the 

matrices 


are those for 

9 

and 

s'" 

put equal to zero- 

If 

the vector 

x<0> 


5"». ana 
except that 0 is 


k( 0) - rv<0> v<0> ir(0> ir<0) >(0^ (A 


then the individual components are given by 


;(0) = (p(l),(l)*^ p(l)%(l))^ , 


(4 


. 21 ) 


. 22 ) 


.23) 


.24) 
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( 0 ) 


= + 0,tw 

L y y 2 


<» - iu<l>*iu‘l>) 


'] 


* (6j + ^ we,)(ip<i>u‘i>* - ip«>*u'i>) 

+ 3y(p'l>v<W* * p(l>*v<l>) 


- [t'{t 


«t/n.(i) _ ie + ie„v*^^) 

y A y ^ 


+ y" Uy 


(4.25) 


( 0 ) _ 






(1)*.-..(1) 


+ 0„(p'-"' iv'-"' - ) 


(1),.(1) 


- {_|7. {T(l>(-i0^u(">* + v^l>* - i02w(">*) 

+ T(l>*(i0,u(l) + + i0,w^">)} 

A y 2 

+ 2 Y,' + T^^^*Vy^^)}y , 


(4.26) 


ti.(0) 

'^4 


y y 


( 1 ) 


+ (0^ + U0JJ + W02> (P'“' - iw 


'»* + ip'i'V") 

+ Wy(p'"v<“* + P<1'V») 

“ |j7' (w^^^* - iOgV^^^*) + + i02V^^^)} 

+ Y" Wy *Jy , (4.27) 
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(0) 


. pc,[e^i(u'l>*T<» - U<»T<1>*) * v<»Tf>* * 




■>*,] 


* (1) <p'»*T<» - p«>t‘»*) 


* p5.f^„(l)^(i)* ^ T<i> V») 






+ 29j((u^» - iv«>* -f »»l*iv'«) 


^ 2v'»v<^>* -P 


P 20.(w'">(-i)v<»* P w»>*iv<«) 

y y 


P (e2p^ej)2»w„«>* P (e^p|6^)2»<i>„«>* 


- §®x®z 


] 


p‘2v«>v<«* - ie,(i»'»v'l>* - iu<“V«) 

3 y y 3 X y y 

- U, Uw<»v<»* - iw<">*v<l')1 
j ^ y y j 


.,(1) 


Uy(Uy^^* - ) + w„(w 




- 2yT^^^*Fu + w + iO.v^^^n 

Lyy ^ yy 


'y"T<1)t(1)*(uJ+wJ) , 


(4 


.28) 
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(4.29) 


_ ( 1 ) {!)* 

Fq2T T 

This problem as defined by (4.22) and the appropriate boundary 
conditions can also be solved at each (X,Z,T') location. 

At this junction f ^ local 

dispersion relation are all kiiown. In the problem for # 

which determines X , all of these functions are required. 

This problem is 

(A, + B, ;^ + C, = 0 , (4.30) 

A/l f\jl dy 0/1 'v 

where 


'U 


,.(11) 

.(11) .(11) ,^(11) fc(ii), 

» ^4 ' *^5 ' ^6 ' ' 

(4.31) 

= 

Or 


,(11) 
^2 ' 

-(11) ,(11) ,(11) ,(11), 

*■3 ' 4 ' 5 ' 6 ' * 

(4.32) 

functions 

3 

and are readily calculated: 


fc(ll) _ 





+ iSj[p 

(2)^(1)* 

+ p<»*w<=> t p<»>w<"> + p<">w"»‘ 

] 



(1) * 
v' ' + 

p(l>*v<2> . pOlpd) ^ p(l)^(0)j^ 

t 


(4.33) 


k<“> = +„W* 2 iu< 2 > - 


+ v< 2 ).jfll* + vl»*„< 2 ) + v«»u'» + v'l’u"" 

y y y y 


+ e^(w«»i«<» - »<2)iu(D* 2i„(l)*„C2),j 


* 5y(p<“'v'" + pWvX” + p<2)v<l)* + p(l>*v<21) 


+ (9,U+ e,S)(ip«»»<l> + 21p<l>*u<2l . lp(2)u<l'*U 




■p e.f-io'i’w'i’u'i’* + ip<i>„<i>*u'i> + ipWVi>u‘i>) 

+ p<l>v'l>u^»* + p‘l>vll>*uW + p(l)*v(l)u<l) 

+ 6j(T<l>*2iw<='> + T«”iw'l> - iT<2)„(l)*) 


- 2 y ' * 2 iu^^^ + 




+.- 2T^^^T*^^*(i0„u^^^ + iO + v*^’)} 

X 2 y 
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- it'” 


+ + 2i0xV^^M} 
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v(ll) 

*3 




+ + w<l>* 


+ »<2>vW* t v«’>v'» + v<“*,<2>} 

y y y 


+ ej{-ip'=>vW - ip<»>v<i> + 2 ip'i'*v< 2 >) 


* e,U(p<»«'i'* p<i'*„'i>)»'*> - ipW„<»v'ii* 


* e,(i(p<‘>w'»* - P<“*w<»)y'*> - ip<^'wWv'l>* 


* tpll>,<l)v^»* + p(l)v'l'*vW + p'l>*v<l'v“') 




'» - ie w'» + v<»>) 

Z y 


+ T<o>(ie,u<"> + i© + v<^>) 
X z y 


+ T(l)*(2i0^u<2) + 2i02w(2) + 


+ 2 [y-(T< 2 )v.</)* + + T<^> V. 2 >) 


y Jy 


']: 


i [Y"{2T^^^T^^^*(i0„u^^^ + i0„w^^^ + 

-j A ii y 


^^(i)^-ivW*-ie,w'»**vW*)) 

. 7., Tilled)* ^ 2T<l)^(l)*,a),-j 

y y I V 


} 

} 
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t fT»>*(w<2> Y 2ie,v<2)) 

y " 

- le,[|.rT'»^«<»* - iv'»*. 

+ + i0 

y 2 

Y i;l 5 y(i)\(i)^ 

2 y J » 




( 4 . 35 ) 
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>(11) _ 
*4 


p[0x(u<^>*2iw<2) _ 

+ 02(iw<2)w(l)* + 

+ e^(2ip<l)V(2) ^ ip(0)„(l) ^ ip(2)„(l)*, 

+ Wy(p^2)^(l)* ^ p(O)^d) ^ pd)^(O) ^ pd)%(2)j 

+ (®x“ '*’ ®z^) (ip 

+ 0x{-ip<">u(l>w(l)* + ip<l>u<")*w(l) + ip(l)*u(l)w<"h 


+ 0,{ip<l>*w(l) } 


+ {p Wy^^ * + P *Wy^^ + p } 

- i0xY‘{0_(iTfO>uflJ - iT<2)u(l>* + 2iT«l)*u<2)j 


+ 0jj(iT + 2iT ) } 


i0jjY^(02 (“iT + 2iT ) 


+ 0„(-iT^^^ + 2iT ) } 


- |y'{T^^^ - ie^v^^^*) + 


+ - 2i0„v^^^ 


)} 


+ Y"Wy{T 
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- + iS.v'W)) 




■f iej[|T'{e^(-lT<2>„<ll* + iTfOlu'!' + 2iT'»*u<^>) 

t e^(-iT<^>w'W* ♦ «'»>»<“ ♦ 2iT<»*w<=>) 




- 27'e^(-iT<2>»<l** + iT<“'w'» -f 2iT<“*„»>)] 


+ ej(-iT<» W<1>* ^ 2iT<l'T<l'*w'l>) 




'•v^n] 

- le^ Y- [-IT<»\'“* - 2tT'“T<l'*w<l>] 


(4.36) 
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+ ej,{-lw<2>T<“* + iw«»T<l> + 2iw'l'*T<2)) 

+ {^(2)j(D* ^ ^(l)*j(2) + (0) (1) ^ v<'->T“'>n 

y y y y J 

+ + “©x ■*■ W02> (2ip 

+ + p(0)v(l> + p(l>v(0) + p(2)^(l)*) 

+ ipc^ (0j + U0JJ + wOg) 

+ Fc^fy (T 

+ + p{l)^(l)*T(l) _ p(l)^(l)^(l)*j 

+ + p(l)„(l)*^(l) _ p(1)„(1)t(1)*) 

+ c^(p * + p *Ty^^ + p ) 

+ i Jc^0jj(2u *T ) + i Fc^02 *T ^ 

+ Fc* *T T + V (T T ) ) 

^ V y ' ^ y ^ 

2 

+ =v^®T ®x“ p^^^* 


+ p *T ) 



(11 * 

PC"(0n. + + 02W)iT^-^^ t'-^' 


I‘’°v"='T 


K. i?S"T (T<W v'»* + 2 t'1>T<1>*v‘») 

z V y 


- y[ 2 u^'’>u 


< 1 > + 

Y y y 


+ 2i0^(u^O>v<l) - + 2 u^1)*v(2) 


u^2)^(l)*j 


+ 2wf>w^l> +.2 w^2)w^1)* 


+ 2i0,(w^v<l> + + 2 w^^>%( 2) . „^2)^(l)*j 


+ 4(0j + 02,v^2)^(l)* ^ ,el^(2)^a)* 


+ 40jw(2)„(D* 


4 - 
3 


- 4 0j^02(u'^'*w'2' + 

7 [40ju(2)u<l>* + 402 w<2)w<1)* + 2(v^2)^U)* 


y y ' 


- ie^(v^“'u'l' - v”>u<»* * 2 v^“*u<">) 


- ie,(v^«>„<» - v^^>w'l>* ^ 2 v^»*w'^>) 


+ 2e^e^(a<2)„<i>* + u'l’V^')] 


27.5-^[t< 2),„»>* - iSj;V<“*i + T<»>(u^i> ^ iejjV<i>) 


* T<»'(«f> t 21 v' 2 )) 


0 
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- 2y*w, 


■[’ 






t(o) („a) 




+T<^>*(w< 2 ) + 2 i 0 v< 2 >) 

y z 


•] 


T-(»^ •f 5j) (T<«>T<» + T<2 't<»*) 




^ 2(»w-ne,v<»HwW* - ie,v<»*) 




+ T<»*((aW * ie,v<i')" + („'1> ^ ie,v<i>)' 

y A y * 


- (e,u'» ^ e^w‘W)')] 

- |v[3jw«>2u'»u<»* - tW*u<»') 

y y y 

+ 0 ^{t' 1 ' 2 w«>wW - T<»*«<»^ 

- lex<T'”»“’v^”* + T<»*uWv<» - T<l>u'W*v<i>) 

- ie,(T<»w<»v<l>* + T<l>*„<»v<i' - tWw<1>*v<1I) 

^ y y y 
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-^(u^ 

+ + 10 

X z y 

+ p‘">*(2ieyu'2) * 2ie^„(2) * ,U),T 

^ <’^02 * f„i)[3''^’(-iv‘”* - i®2>'‘^’* + V«>*) 

+ T^®^(i0 + 10 

A z y 

+ + 2102W^2) ^ 

+ T{|f2iP^^^ + 

{-10„u(D* - 10 

A z y 

+ T{F2iP<^^p(^>* + + p(l)*T^^>) + Fq3t(1)t(1>*} 

{lOjjU^^^ + + Vy^^} 
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+ {Pll(p<^>T<l>* + p<1>*t(^>) + 2Fo2T<1Jt<^>*} 


{ i0„u^^^ + ie„w^^^ + 

X Z y 


+ k*(eE + 


X z 




yy 


|k"(02 + 02)T.<") 


- |(k-T<^> , 


(4.37) 




F^^(P< 0 )t< 1 > + p(l>TtO> + p< 2 )t< 1 >* + p(l>*T( 2 ), 

Fo2<T^°^T^^^ + 

|f3qP^^^ P^^^* - 1,^21 + 2p^^^p^^^*T^^^) 


ip^,W'»%'l>* * 2T<»T<»*p'l') - ip„3T<l'\<»* 


(4.38) 
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J 



,( 11 ) _ 
h 

,( 11 ) _ 






( 11 ) - 


( 11 ) 






( 11 ) _ 


, (4.39) 


£^ 11 ) = 0 


(4.40) 


.For the nonparallel contributions we must examine the 
series of problems 

dy^ 


(^1 TT ^1 = 0 ' 


j = 1,2,3, 


(4.41) 


where 


K<1> 

- 

' ji * 32 ' j3 ' j4 ' j5 ' 

’"j6 ^ ' 

(4.42) 


_ ,.(1) ,(1) .(1) ,(1) »(1) 

*"6 ^ * 

(4.43) 

These functions r twenty-four of them in all. 

are now listed. 


^11 

= p(l>u + PU<1> , 


(4.44) 

^12 


+ + i0_w^^^) 

y ^ 



- 27 i 0 jjU^^^ + - 

2Yi0x'^^^^ 



- (Yv^^^)y - , 

• 

(4.45) 
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^13 


- yQ + 10 

A y A 


- Y'U + t (YU^^^)„ , 


(4.46) 


^14 




- Yi©x''^^^ ■•■ 


(4.47) 


= 

'l5 


pc^uT^^^ - 2YUyV^^^ + 


- ik©jjT^^^ - IGjjkT^^^ , 


(4.48) 


^(1) 

^16 


0 , 


(4.49) 


( 1 ) _ 




( 1 ) _ 


pu<^> 


( 1 ) _ 


pv^^^ , (4.50) 


(1) - 




( 1 ) _ 


PC,T<1) 


( 1 ) _ 


= 0 


(4.51) 


c<l) = 
^21 


= P 


+ pw(l> , 


(4.52) 


= pwu^^^ + l-Ye..iw^^^ - Y0_iu^^^ 


'22 


■^ye^iw - Y^g^*^ 


- + 0jjW^^^) , 


(4.53) 
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I 


^23 




- + i0„v^^^) - y’w„T^^^ , 


(4.54) 


= pww^^^ - + P^^^ 


"24 


+ |y[10xU 


+ i0 


i0,w<^>] 


- 2Y02iw^^^ - (YV^^^) 


+ “ 2Y0^iw^^^ , 


(4.55) 


5^25^ = - 2YWyV^^^ + 


- k02iT^^^ - 02kiT^^^ 


(4.56) 


= 0 • 


(4.57) 


.(1) = /- 
"31 


= («x ■*■ "^y '*' 




+ i(u0jj + W02)p^^^ + ip(0jjU^^^ + 0gW^^^) 




(1);t ^ T..(l) 


+ (PV' My + Px “ + P*^X 


^ p»'i; . fw<» , 


(4.58) 
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I 


+ + uujj + vUy + wu2)p^^^ 

+ jY'TjjUejjU^^^ + + iOgW^^^) 

- 2Y0xxi’i^^^ - 2Y’Tjji0^u^^^ 

t |7'V’'‘"<“x * ^y * "z> 

- 2Y'0xUxi'’^^^^ ~ iY’OgT^^^ (u^2 + w^^) 

- Y(©zziu^^^ + 

^ Y'TgdGgU^^^ + 

+ i{(pu + pu) 0 jj + (pw + pw)02>u^^^ 

+ (pUy + pUy)v^^^ 

- |•Y'T0x(0x'^^^^ + + 2 y'T0^u^^^ 

+ y'T 02(02U^^^ + ©xW^^^ - (y'TCu^^^ + i0jjV^^^}y 

- ((Y'Uy + Y"UyT)T^^^)y +• puuj[^^ + p^^^ 
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k 


( 1 ) 

33 


+ |Y(i0xU^^^ + 


10 ,.-o „(1) _ ..(I) 


,n-„(l) 


- ^ iY0xU' ' - iY02;W' ' - (yv'"'). 


+ PVTU^^^ - |iY0xW^^^ , 

pVyV^^^ - Y0xx^^^^^ “ Y'Tjj(Uy^^ + 

- (Y'Uy)jjT^^^ + I ( Y* (Ujj + Vy + W2)T^^My 

- 2(7'VyT^^^)y - Y'TgCw^^’ + 102 

- - (7'Wy)2T^^^ 


+ i{(pu + pu)0jj + {pw + 

pw)02>v^^^ 

- 1y'T(u^^ + 10 

y ^ 


- 10jj(Y'Uy + Y”Uy?)T^^^ 

- 10„(y'w„ + Y"w„T)T 
^ y y 

- 1y'T02(w^^^ + 102V^^^) 

- 2(Y’Tv^^^)y 

- |^y’T(10„u^^^ + 

^ /i. y 

+ 102W^^^) 


+ puv^l> - i9xYvf> - T(«^” + ivi“> - 


(4 


(1) 


’fVx" 


‘Vr-a. 


.59) 





- |(Tu(^^)y - 

" y^'^yV ^ , (4,60) 

1^34^ = + pWgW^^^ + + uwjj + Wy + wwg) 

- Y’TjjdOgU^^^ + i0jjW^^^ - Y(02xiu^^^ + 

- Y’iQx'^^^Nug + Wjj) +1 Y'T2(i0jjU^^^ + + i02W^^^) 

+ +i0Zz"‘^^> - 2(Y0z)ziw<^> 

+ |Y0ziT^^^ (Ujj + Vy + Wg) - 2 y'02Wz^’’^^^ 

+ { (pu + pu)0„ + (pw + pw)0_}w^^^ + (pw„ + pw„)v^^^ 

A /» y y 

- |y'T02(0x“^^^ + 0z”^^^ “ 


{y’T(w^^^ 

+ i02V^^^}y 



2y'T02W^^^ 

- ((Y'Wy + 

Y"TWy)T^^^)y 


{y'T(w^^^ 

+ i0zV^^^}y 

+ y'T0jj(0jjU^^^ 

+ 02« 

■pUW^^^ :- 

ir<iSz"iS" i®x"x“> - wrf’ 



I 


+ - 4i0„ywi^^ 


rz 


'Z^«Z 


. |r[io,uf> * v^. . ie,w<»] - 


(Yv‘») 


yz 
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k||^ = pc„{yT^,^^ + 


V y 


y[2i(u2 + Wjj) 


+ j + 2VyV^^^ + 2iw202W^^^ 


<"x''y" * i®xV‘”> 


(WjV^l> + iejVy«<^>) 


+ T (Ujj + Vy + w^) + ^ 02 '“^^^^^ 


* ^01 ■^7y 


- ik* (Tj^0jj + T202>T^^^ - i{(k0jj)jj + (k02)2>T^^^ 


,-n,(l) 


^■^\ A rs rp(l) 


+ (pc^, + pc^T)0^iT' + (pc^u + pc^Tu + pc^u)i0yT 


V " 


+ (PC^Ty + PC^TTy + pc^Ty)v 


(1) 


- (2yUy + 2y'TUy) ) 


- (2yWy + 2y'TWy) + i02V^^^ 


2r‘(SyUy + WyWy)T 


_ Y"T(Uy + Wy)T^^^ 


Y’(Uy«y + WyWy)T 


( 1 ) 


- Y"T(Uy + Wy)T^^^ 


+ (TPqj^ + T(Pj^^p + Fq 2T)) 


+ ic'T {0^ + 

A. li 


(k'TT^^^)y 


rz Zn,(l) _ orz „(1) 


( 1 ) 


+ pc^uT^-"' - 2 yu„v^"'' - - 2ik0^T 


.(1) 


y X 


D1 X 


X X 




(4.62) 


,(1) = 
^36 


<^20P 


{F^^p + Fq2T)T^^^ 


(4.63) 


The double bar symbols in these equations correspond to 

the first correction to the basic boundary layer flow* This com- 
pletes the detailed listings of these functions. 

It remains to show how X and j = 1,2,3 are to be 

calculated. For this purpose the return to the basic linear 
eigenvalue problem associated with equations (4 . 3 , 4 , 5 , 6, 7 , 8) . 

This can be written as 


dV 

Lv = Drt — ^ 3 ^ + D V = 0 

'X/ 'vO ^^2 1 dy 'v2'\^ 


(4.64) 


subject to the appropriate boundary conditions. Here v is the 


vector (p 


(1) „(1) „(1) ,(1) ^(1) ^(1) 


u 


, w'^^ , T' 


p ' ' ) and D. j = 0,1,2 
A/3 


are 6x6 matrices. 


52 



The adjoint problem is that associated with the equation 


^ ‘SoS> - * & - 0 ■ 


(4.65) 


and the adjoint boundary conditions follow the identity 


T T 

(w Lv - V Mw)dy 

f\j <Vi 'Vi 


fp dv m 'T — n OO 

" ^ (4.66) 

0 dy A/ 'bl dy 'x, '\,0 


It follows that the inhomogeneous problem 


Lv = f 

Or 


(4.67) 


will have a solution if and only if 

rOO 

w^fdy = 0 . 

Jo 


(4.68) 


On applying this condition to equations (4.30) and (4.41) 
we obtain 


X = - 


w^K<^">dy 

Or ^ ^ 


w’'L<ll>dy 

0 -V 'v. 


(4.69) 


w’^L<">dy 


j = 1,2,3 


(4.70) 
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with these functions determined at each location the 
coefficients in equation governing the amplitude A 
(equation 4.12) are known and the spatial and temporal growth 
of A can be studied. 



5. CHEBYSHEV- SPECTRAL METHODS POR STABILITY CALCULATIONS 


In Sects. 5-8, we discuss the state-of-the-art of 
numerical techniques for staibility calculations. In this 
Section, we provide an introduction to the use of Chebyshev 

9 

polynomials for the solution of stability problems. 

An important difference between finite-difference 
approximations to the eigenfunctions and eigenvalues of 
a stability problem cind Chebyshev polynomial approximations 
to the same problem is their order of accuracy. Finite- 
difference approximations give only a finite order of accuracy 
in the sense that errors behave asymptotically like h^ for 
some finite p when the grid spacing h approaches zero. On 
the other hand, if the unperturbed velocity profile is smooth 
(infinitely differentiable) , the Chebyshev polynomial 
approximations discussed here are of infinite order in the 
sense that errors decrease more rapidly than any power of 
1/N as N » . 

Another difference between finite-difference and 

Chebyshev polynomial approximations to stability problems 

concerns their resolution of possible regions of rapid 

change {'boundary layers*) in the eigenfunctions. When 

the Reynolds number R (based on boundary layer thickness 6 

and frees tream velocity U) is large, the eigenfunctions 

- 1/2 

exhibit boundary layers of thickness of order R ' near 

-1/3 

y = 0 and internal layers of thickness of order R ' near 

the critical layer (where wave phase speed equals flow velocity) . 



In order for finite-difference approximations to be 

accurate, it is necessary that the grid spacing be at most 

near y = 0 and at most near the critical layer. 

Thus, if uniform grid spacing is used the number of grid points 

1/2 

required is scaled by a factor R as the Reynolds number 

increases. If non-uniform grid spacing is used, this difficulty 

may be partially avoided. On the other hand, the number of 

Chebyshev polynomials required^ for accurate stability 

1/4 

calculations scales only as R as the Reynolds number R 
approaches infinity. In laminar flow control applications, 
this difference between finite-difference eind Chebyshev 
polynomial methods is important. If the range of Reynolds 
nxmbers to be studied in a given LFC application is, say, 

R = 1000 —10 ,000 .and if, say, 20 polynomials or 100 grid 
points are required to solve the problem at R =1000 (these 
resolutions are, in fact, typical), then less than 40 polynomials 
will be required at R=10,000 while more than 300 grid points 
will be required. 

The rapid convergence properties of Chebyshev polynomials 
are verified as follows. If the unperturbed flow is smooth, 
then so are the eigenfunctions of the linearized Navier-Stokes 
equations. Let T^ (x) denote the nth-degree Chebyshev polynomial 
of the first kind, defined by 

Tj^(cos 0) = cos n0 (5.1) 

for all non-negative integers n. Some examples are (x) = 1, 

• 2 

T^(x) = X, T2<x) = 2x -1. It is possible to expand the 
eigenfunction ^(y) in the interval -1 £ y £ 1 (we discuss in 



Sec* 6 techniques for handling the semi- in finite interval 


0 y ^ encountered in boundary- layer stability problems) as 


4^(y) = 2 (5.2) 

n=0 ^ ^ 


where 

= il“ I ’f'{y)T„(y) dV) dy (5.3) 

n ^ —I 

with Cq = 2, ^ n >0. The rapidity of convergence 

of (5.2) for |y| ^1 is easily demonstrated by observing that 


f(0) = iMcos 0) 

is an infinitely differentiable^ even, periodic function 
of 0 . Consequently, the theory of Fourier series ensures 
that f(0) possesses a Fourier cosine expansion 


00 

" n£o ^n (5.4) 

with the property that the error after N terms decreases 
more rapidly than any power of 1/N as N The expansion 

(5.4) is precisely (5.1) for y = cos 0 . 

The infinite-order accuracy of Chebyshev polynomial 

approximations to smooth fiinctions holds no matter what 

the boundary values of the functions or their derivatives, 

in contrast to the situation when other classes of orthogonal 
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functions (like trigonometric or Bessel functions) are used. 

In the following stabsections, we discuss several programming 
and technical aspects of the application of Chebyshev polynomials 
to stcibility calculations. 



Chebyshev Matrix Method 


The derivation of the equations satisfied by the 
expansion coefficients in the Chebyshev expansion 

N 

’My) * 2 a T (y) 

n=0 " " 

can be difficult and time-consuming if' done by hand. In 
order to improve the flexibility of this method and the 
ease in which it can be applied to new problems, we have 
developed a nearly automatic matrix method for computer- 
generation of the equations satisfied by the coefficients 


To illustrate the method, suppose we defined the 
vector jp by 


If, = 




/ 


(5.5) 


(5.6) 


When the expansion (5.5) is substituted into a differential 
operator, equations for the coefficients a^ are obtained by 
re-expanding the result in series of Chebyshev polynomials 
and equating coefficients of each Chebyshev polynomial: if 

- nloVn'y> 

then each b is a function of the a and we obtain N equations 
n n 

approximating dC ^ = 0 t>y setting = 0 for n £ N. In order 
to find these equations we must derive an efficient and easy 
method to determine the effect of the differential operator cjC 
on the vector 



Consider the Chebyshev representation of the function 


♦•(y) - ^ 


Let the associated vector of Chebyshev coefficients of i|;* (y) 


be denoted by 
Since 

+ V3 + •••>' 

it follows that 




( 1 ) = 2 


pa. 


" ‘"n p=n+l 

p+n odd 

where Cq = 2, = 1 for n >0. Therefore, 

I* = 

where the N+1 x N+1 matrix D has elements 

= 0 ifi>jor i+j is even 

otherwise 

°i-l 

Similarly, if f (y) = y\p (y) , then the Chebyshev coeffi- 
cients of f are' given by 

-> 

f = Yip 

where the N+1 x N+1 matrix Y has elements 
Y^j = 1 if i=2, j=l 

1/2 if i=j+l or i=j-l 
0 otherwise. 


The utility of setting up these matrices (and other 
similar ones) is the ease with which they may be used to 
set up the matrices expressing complicated differential 
operators. For example, the matrix for the operator 



yd^/dy^(y^d/dy) is just YD^Y^D. (In fact, the last 

2 3 

statement is not quite true because the matrxx YD Y D 
does not correctly represent the action of the differential 
operator on high-order Chebyshev polynomials. However, it 
may be shown that this error is nearly negligible and we 
will not discuss it further here.) 

In actual practice, it is frequently not necessary to 
store the matrices D, Y, etc., because their veiy simple form 
makes it possible to generate their elements as they are 
needed during the computation. 

Let us consider how the Chebyshev matrix method applies 
to the solution of the Orr-Sommerfeld equation 

(dV<3y^ - a^)^l'(y) = iaR(5(y)-c) (dVdy^ - a^)<My) 

- icHu” (y)'l’(y) , 

where u(y) is the unperturbed profile. In terms of 
Chebyshev matrices, these equations are 

(D^ - = iaR(U - cl) (D^ - a^I)| - iaRU”| (5.7) 

where U is the matrix that multiplies by u(y) and U" is 
the matrix that multiplies by U” . 

Tau Method 

The equations (5.7) for the Chebyshev coefficients 
a^ do not account for the boundary conditions imposed 
on iMy) . There are several ways to impose the boundary 
conditions consistently on (5.7)? this can be done by 
Galerkin, collocation, or tau approximation^’^ It is 
usually most convenient to apply tau approximation, as 
we will now discuss. 



The idea of the tau method is to drop enough of 
the equations (5.7) that all the boundary conditions can 
be applied. For the Or r- Sommer f eld equation, the boundary 
conditions 

iKy) = (y) = 0 

should be applied at the rigid boundaries y = ±1. Thus, 
we delete the last four rows of the matrix equation (5.7) 
and replace them by the four boiandary conditions 


N 

E 

n=0 


n 


N 

2 

n=0 


(-1)" a 


n 


0 


No N _ , -1 o 

Z n a„ = Z (-1) n = 0 

n=0 " n=0 ^ 


(5.8) 


Eqs. (5.8) follow because T^(±l) =(±1)^, T^(±l) = n^(±l)^^^. 
Thus, we retain N-3 equations of the form (5.7) and 4 equations 
of the form (5.8) so that there are N+1 equations for the 
N+1 \mknowns aQ,a^,...,a^. 


Fast -Fourier Transform 

An additional advantage of Chebyshev polynomials over 
many other orthogonal bases is the existence of the fast Fourier 
transform to effect efficient conversion between the Chebyshev 
coordinates a^ and the physical space perturbation ^(y). In 
fact, since 

N 

(cos 0) = 2 a cos n0 , 

n=0 " 

Chebyshev series can be summed by any technique that sums 
Fourier cosine series, in particular, certain variants of 
the fast Fourier transform (see the Appendix of Ref. 12). 



6. NUMERICAL METHODS FOR EIGENVALUE CAIiCULATIONS 


We begin by making a number of general comments that 
apply to most niimerical schemes for stability calculations. 

First, we remind the reader that the problem is difficult 
because the Or r- Sommer f eld equation (or other linearized equation) 
is moderately stiff at large Reynolds numbers, as reflected by 
the boundary layers exhibited by eigenfunctions (see Sec. 5) 
Second, calculations of moderately high accuracy are 
frequently needed for such purposes as computation of the 
group velocity (see Sec. 7) and various optimization schemes 
(see Sec. 7). Third, we comment that there are two possible 
kinds of stability calculations that can be made, spatial 
stcibility and temporal stability (see Sects. 2-3). For 

small growth rates, the results of temporal and spatial stability 

13 

analyses are closely related^ However, because the differences 
between these two types of analyses may be central to the 
problem of LFC aircraft design, let us contrast them briefly 
here. 

In a temporal stability analysis, the wavenumbers of 
the disturbance are assumed real while the frequency of the 
disturbance may be complex (a positive imaginary part indicates 
instability) . On the other hand, in a spatial stability 
analysis, the frequency .is taken real, while a positive 
imaginary part in a wavenumber indicates instability. 

Since many problems of aerodynamical interest are, on average, 
stationary in a suitable coordinate frame, it seems that 



spatial stability analysis should be the more relevant one. 

However, this is by no means clear experimentally. It seems 

that low frequency disturbances are, treated better by spatial 

theory than by temporal, but high frequency disturbances 

14 

seem to agree better with temporal theory. 

The confusion between spatial and temporal theory 
is even more severe in the case of the propagation of 
three-dimensional disturbances, which are of primary interest 
in the boundary layers of LFC aircraft. Spatial stability 
theory is ambiguous for three-dimensional disturbances. 

In the case of two-dimensional disturbances, it is physically 
plausible that the direction of maximum growth of the 
disturbance is perpendicular to the constant-.phase surfaces 
of the disturbance and parallel to the freestream flow 
direction. On the other hand, in three-dimensional layered 
flows, there is no apparent reason why the direction of 
maximum growth should be perpendicular to the direction 
of the constant-phase surfaces. If these directions are 
allowed to be arbitrary, one quickly gets involved in 
ill-posed mathematical problems. Until this basic question 
is resolved, it may be best to use temporal stability theory 
and a group velocity transformation for three-dimensional 
stability analyses of LFC boundary layers. 

Finally, we comment on the mathematical technique to 
treat the semi-infinite domain 0 <_ y < » of the boundary 
layer. Grosch and Orszag^^ have shown that the best way 
to handle the y-direction is to transform it by means of 
the algebraic transformation 


into the finite interval -1 ± Y < 1, and then apply standard 
numerical techniques on the transformed finite interval. 

This technique provides accurate results with about 50% 
less resolution than required using simple truncation of the 
domain 0 £ y £ L. Some results of a Cheby she v- spectral 
calculation are given in Table 1. Here we compare the 
accuracy of the most unstable mode of Blasius flow at a 
Reynolds number of 580 and a wavenumber of .179 using several 
treatments of the boundary at y = » . The methods are: 

(i) truncation, which involves solving the problem on the 

finite interval 0 £ y £ L for several values of L and with 
both no-slip (L) = (L) = 0 boundary conditions and 

asymptotic boundary conditions tJ/* (L) + aifiCL) = 0 applied; 

(ii) an exponential map of the form Y = 1 - 2 exp(-y/L); 
and (iii)the algebraic transformation (6.1). For both 
methods (ii) and (iii) , several kinds of boundary 
conditions are applied at Y = 1 (y = " ) , including no-slip 
conditions and boundary conditions ( ! ) . We conclude from 
Table 1 that no boundary conditions at all need be applied 
at Y = 1. 

With the mapping (6.1), L may be chosen to optimize 
the accuracy of the calculations. After some experience, 
it has been found that a good choice for h ±s h ^ ^^0' 
where y^ is the value of y at which the streamwise component 
of the velocity achieves 1/2 of its freestream value. 



Table 4, Eigenvalues of the Orr-Sommerfeld equation for Blasius flow, R = 580, a « 0.179 


Case 1 

Mapping 

Boundary Conditions 
at 'z = 

L 

N 

(Number of 
Chebyshev Modes) 

c 

i 

!■ 

Truncation 

4»(L) = (L) = 0 

10 

1 

44 

0.37887 7 + iO. 00025 0 1 

2 


II 

20 1 

44 

0.36455 7 + iO. 00777 3 i 

.3 


II 

20 

1 

46 

0.36455 1 + iO. 00778 1 = 

4 


« 

; 

30 

44 

0.36399 6 + iO. 00788 8 | 

5 


(L) + at{»(L) = 0 

20 

44 

0.36021 3 + iO. 00667 1 [ 

6 


II 

30 

44 

0.36404 1 + iO. 00811 3 ! 

7 

Exponential 

’l'(l) = 0(1) 

1 

42 

0.34858 0 + iO. 01312 9 

8 


II 

1 

46 

0.34961 1 + iO. 01285 6 

9 


I|»(l) = ( 1 ) = 0 

1 

46 

0.38378 9 - iO. 00276 6 

10 


II 

1 

1 

70 

0.37853 1 + iO. 00047 1 

11 

Algebraic 

<|»(1) = 0(1) 

1 

26 

0.36414 7 + iO. 00800 7 

12 


II 

1 

34 

0.36412 1 + iO. 00795 76 

13 


II 

1 

42 

0. 36412- 288 + iO. 00795 975 

14 


4»{1) = (1) = 0 

1 

42 

0.'36412 325 + iO. 00785 894 

15 


II 

1 

60 

0.36412 285 + iO. 00795 973 

16 


<|)(JL) =0 

1 

42 

0.36412 287 + iO. 00795 976 


c\ 

ui 



Global Methods for Temporal .Eigenvalue Calculations 


When no guess is available for the eigenvalue of 
interest r it is best to use a method that is globally 
convergent and nearly guaranteed to converge to the 
eigenvalue. Such a method may be based on the matrix QR 
algorithm^^ for calculation of the eigenvalues of a 
general complex matrix. 

When the Orr-Sommerfeld equation (or other similar 
differential equation) is formulated as a matrix problem 
(using either Chebyshev polynomials or finite-difference 
methods ) , it takes the form 


^ =? ABx 


( 6 . 2 ) 


where A is the eigenvalue (denoted c or m above in the 
case of temporal stability calculations) eind x is the 
discrete representation of .the eigenfunction. The eigenvalue 
is determined by the determinant condition 

Det (a - XB| = 0. (6.3) 

Eq. (6.3) is a generalized eigenvalue problem and the 
matrix QR algorithm does not apply directly unless either 
A or B is invertible. (If^ say, B*^ exists then (6.3) is 
equivalent to the standard eigenvalue problem 
Det |b“^A - l| = 0.) 

However, it is frequently the case that A and B are singular 
and a more general method must be developed. 

To solve the generalized eigenvalue problem with singular 
A and B, we proceed as follows. There are two steps that are 



executed recursively. 

(i) We use fully pivoted row operations to reduce B to 

upper triangular form, executing the same row operations on 
the matrix A. The resulting generalized eigenvalue problem 
Det I A* - 0 has the same eigenvalues as (6.3) because 

the same row operations were performed on A and B. If airithe 
diagonal elements of B* are nonzero then B* (and hence B) is 
nonsingular and the problem Ccin be immediately reduced to the 
standard eigenvalue problem (and then solved by the QR algorithm) • 
Thus, let us assume that all elements b|j with K £ i ^ N are 
zero (if any elements of this matrix were not zero, full pivoting 
would ensure additional nonzero diagonal elements) . 

(ii) We perform fully pivoted column operations on the rows 
j=K,...,N of the matrix A* to transform A* into an upper 
triangular matrix A" . The same column operations are performed 
on B" but it is still true that b?j = 0 for K £ i £ N and all j. 
The generalized eigenvalue problem Det |A" - XB" ] = 0 has the 
same eigenvalues as (6.3) because it is obtained from it by 

row and column operations simulatneously on both A and B. 

However, rows K,..., N of A" - XB" are upper triangular, 
so the generalized eigenvalue problem for A" and B" has a 
solution only when the generalized eigenvalue problem 

Det |a^ - XB^I = 0 

has a solution, where A^ and B^ are the K-1 x K-1 dimensional 
matrices obtained from A" and B" by discarding rows and columns 
K, K+1,..., N. 

(iii) Go to step (i) . Eventually B* must be non-singular or 
there are no generalized eigenvalues or all 1 are generalized 
eigenvalues • 
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Some results of using the above algorithm on a CDC 7600 

•f* 

computer are given in Ted^le 2. ‘ 

Table 2> 

Timings of the Global Eigenvalue Program 


26 X 26 43 X 43 

matrix matrix 


QR method 0.15 s 0.55 s 

(Fortran, 
all eigenvalues) 

Generalised 

eigenvalue problem 0*18s 0.66 s 

(invertible A or B) 

Generalized 
eigenvalue problem 
(singular A and.B) 

Two recursions 0.27s 1.1s 


The advantage of this technique is that it is 
very general and is globally convergent. The disadvantage 
is that it is not too fast (though it is probably faster 
when combined with Chebyshev polynomials than most commonly 
used methods) • 

When finite-difference methods are used, it may be 
better to obtain a globally convergent result by use of 
the LR algorithm and band matrix factorizations. We have 
not yet tried this technique and cannot yet comment on it. 

much simpler technique to solve the generalized 
eigenvalue problem has been suggested by J. Shearer 
(private communication, 1977) • If y is not one of 
the eigenvalues of (6.2) then A - yB is invertible. 

If the eigenvalues of (A-yB)"lB are denoted by c, then 
the eigenvalues of (6.2) are X = y + 1/c. 



Spurious Unstable Modes 


One of the drawbacks of the global method as formulated above 
is that the generalized eigenvalue problem (6*3) may indicate 
the existence of growing (unstable) modes that are not 
physically relevant. These spurious unstable modes, which 
may appear for either the Chebyshev-spectral or finite- 
difference methods, do not correspond to solutions of the 
differential equation — as the spatial resolution used to 
discretize the eigenfunction changes (i.e. the number of 
grid points or Chebyshev polynomials) , true modes of the 
differential equation converge while spurious modes do not. 

A clumsy way to distinguish spurious modes from true 
modes is to change the spatial resolution and retain only 
those modes that do not change appreciably. This is neither 
efficient nor elegant. 

A better way is to eliminate the spurious unstable 
modes entirely. Spurious stable modes are still possible, 
but since these stable modes are normally very stable, they 
are not of much interest and can be easily disregarded without 
testing their true nature. We shall now describe a technique 
for eliminating the spurious unstable modes. 

The idea is simply that the spurious unstable modes 
would, if we used the same numerical method used for. 
the stability problem on an initial-value problem instead, 
cause the unconditional instability of the numerical solution 
of the initial-value problem. On. the. other hand, if we were 
careful enough to use a numerical method for the stability 
problem that was also niimerically stable for the initial-value 



problem, then no spurious unstable modes would exist. 

There are several ways to eliminate the spurious 
unstable modes in the Chebyshev spectral methods outlined 
in Sec. 5. One way is to use a Galerkin procedure instead 
of the tau procedure discussed in Sec. 5^"^ Another way 
is to factor the fourth-order Or r- Sommer f eld equation into 
two second-order equations and then apply the tau method 
to each of the second-order equations. Thus, the usual 
procedure is to apply four boundary conditions to the 
fourth-order Orr- Sommer f eld equation 

This procedure usually leads to spurious unstable modes.. 
However^ rewriting the Orr- Sommer f eld equation as the 
two second-order equations 

5 = 

" Vy = 


2 ind then applying the boundary conditions = 0 to 
the first equation (because this first equation is equivalent 
to the incompressibility constraint and the associated boundary 
condition should be zero normal flow) and = 0 to the second 
equation (because the second equation embodies the viscous , 
frictional effects and the associated boundary condition is 
no-slip) gives no spurious unstable modes. 



Local Methods for Temporal Eigenvalue Calculations 


There are many local methods (in which a reasonable 

guess is available) for eigenvalue problems* Among the 

well-known methods that have been implemented for difference 
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methods are the methods of orthogonal! zation and parallel 
18 

shooting together with either a bisection search or a 
Newton's method search* There are several effective 
computer codes that implement these procedures, including 
the SUPORT code (written by Scott and Watts of Sandia Laboratory) , 
the TAPS code (written by Gentry and Wazzan of McDonnell Douglas 
Corp.) and codes by Mack of Jet Propulsion Laboratory and 
Keller and Cebeci. At the end of this section, we will 
cite some experience we have had with the SUPORT code cind 
quote some private communications concerning the efficiency 
of the other codes. 


Another way to perform a local analysis is to use 
a simple iterative method to find the eigenvalues of the 
matrix equation (6.3) that approximates the Orr-Sommerfeld 
equation. An effective and efficient procedure for doing 
this is to use the inverse iteration procedure^^ : 

(A - = cx,^ (6.4) 

(A - = c'Yj^ (6.5) 

^k+1 = 


The procedure (6.4-6) is very effective once a good guess 
for an eigenvalue is available because the convergence is 
cubic: 


Here c and c* 


Ajj^^ - A = 0((Aj^ - A)^) . 

* 

are normalization constants so Xj^ 


and yj^ are normalized. 


The method (6.4-6) is designed to work efficiently even 
with non-syinmetric matrices A. Eq. (6.4) shoudl be solved 
using a fully pivoted LU method^ in which case the same LU 
factorization applies to the solution of (6.5). 

In practice, it is not necessary to update the eigenvalue 
approximation after each iteration. In fact, we have 

found it to be most efficient to iterate (6.4-5) approximately 
5-10 times while keeping fixed (and, therefore, using 

the same LU factorization of A - • 

The generalization of (6.4-6) to the generalized 
eigenvalue problem (6.3) is: 


(A - 


= cBx,^ 

(6.7) 

(A - 

^k®^ ^^k+1 

T 

- <= ® s-k 

(6.8) 


1) 

H 

'“k''''^*k 

(6.9) 


The algorithm (6.7-9) still has a rapid rate of convergence 
because it is equivalent to Newton’s method for the solution 
of the generalized eigenvalue problem (6.3) (see below). 

Some data on the speed of the algorithm (6.7-9) for Chebyshev 
methods applied to the stability of interior flows is given 
in Table 3. 


Table 3. Timings of the Local Eigenvalue Program 


26 X 26 43 X 43 

matrix matrix 


Fortran program 0.03 s 0.13 s 

1 eigenvalue 

good guess available 

final accuracy 10““® 

Assembly LU program 0.02s 0.08s 

otherwise same 
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Spatial Stability Calculations 


The Orr-Sonmierfeld equation for spatial stability 
calculations involves a nonlinear, quartic polynomial 
eigenvalue problem of the form 

A(X)x = (X^A^ + X^A^ + X^A 2 + X^A^ + Aq)x = 0 (6.10) 

Global methods for the solution of nonlinear eigenvalue 
problems like (6.10) may be inefficient. A simple global 
method is to set 


= X 

X2 = Xx^ 

X3 = XX2 
X4 = XX3 


A3X4 + A2X3 + A^x^ + AqX^ 


4 4 


and then formulate (6. 11-14) as a generalized eigenvalue 
problem of the form (6.3) involving 4N x 4N matrices: 




(6.11) 

( 6 . 12 ) 

(6.13) 

(6.14) 


The obvious disadveintage of this method is that it requires 
16 times as much memory as (6.3) requires and 64 times as much 
coir^uter times. It is not practical except for problems of 
extreme urgency, in which case it has the virtue of easy 
programming (so long as N is not too big) . 


A local method may be based on Newton’s method for 
the solution of the nonlinear coupled system of equations 


A(X) X = 0 
A(A)’^y = 0 

y^x = 1 


(6.15) 

(6.16) 

(6.17) 


so that y is the adjoint eigenvector to x. Eqs, (6.15-17) 
represent 2N+1 equations for the 2N+1 unknowns , 

y^,X . Newton’s method involves linearizing (6.15)- (6.17) 
about an approximate solution 

A(X)x « A(X^°^)x^°^ + A(X^°^ (x-x^°^) + A* (X^°^x^°^ (X-X^°^ 


= 0 


(6.15*) 


yTx = y(0)T^(0) ^ y(0)T(^_^(0), ^ (y.y(0))T^(0)^ 

+ . . • = 1 

Thus, if we insist that the new approximation x^^^ , X 

be such that the linear approximations (6.15*-17*) are satisfied 
then we obtain the following iteration scheme: 


A(Xj^)Xk+l - 

CA* (Xk)Xk 

(6.18) 

A(Xk)’’yk+i = 


(6.19) 

^k+1 ^ \ ~ 

^k+1 *'*^k*^k 

( 6 . 20 ) 


The advantages of this method are that it is essentially 
as fast as the local method for the linear generalized eigenvalue 
problem (6.7-9) and that it has low memory requirements. The 
disadvantage is that the initial guess sometimes has to be quite 
good for it to work. 

Our recommendation for the most efficient technique 
for spatial stability analysis is as follows: 



(i) If no good approximation to the eigenvalue is available^ 

perform a temporal stcdDility analysis using a globally 
convergent algorithm* Then calculate the group velocity 
using the methods to be discussed in Sec* 7 and transform 
this temporal mode into an approximate spatial mode using 
Caster's group velocity treuis formation^ ^* 

(ii) If a good guess is available (say by method (i)) then 
use the local algorithm (6*18-20) to improve the 
eigenvalue approximation. 


Comparison of Numerical Methods for Plane Poiseuilie Flow 

In this subsection y we present some numerical results 
concerning the stability of plane Poiseuilie flow* In Table 
4, we list the Chebyshev approximation to the most unstable 
mode of plane Poiseuilie flow at R = 10,000 with wavenxxmber 
a = 1 as a function of the number of retained Chebyshev 
polynomials (here we have assumed that the mode feeing sought 
is symmetric in y so we actually retain Chebyshev polynomials 
up to degree 2M) * It is apparent that accurate results 
are achieved rapidly as the number of retained polynomials 
increases* Similar results have been obtained at Reynolds 
numbers of 10^ and higher; at R = 10^, 50 Chebyshev polynomials 

4 

yield *an eigenvalue accurate to about 1 part in 10 (this 
calculation using a global code requires less than 1*5 s of 
CDC 7600 time) * 

In Table 5, we list some results and computer times 
obtained using the SUPORT code mentioned above together 
with a Newton's method iteration procedure for the eigenvalue* 
Although we do not have the data, similar results and computer 
timings have been reported by other groups using finite- 



Table 4 . Chebyshev approximation to the most unstable 
mode of plane Poiseuille flow for a = 1, R = 10000. 

M + 1 X 


17 

.23713751 

+ 

£.00563644 

20 

.23752676 

T 

£.00373427 

23 

.23752670 

+ 

£.00373982 

26 

.23752648 

+ 

£.00373967 

29 

.23752649 

+ 

£.00373967 

38 

.23752649 

+ 

£.00373967 


50 


23752549 + £.00373967 



difference codes. Mack reports that it requires about 3 s 
of UNIVAC 1108 time to converge to an accurate eigenvalue 
given a reasonably good guess; Keller and Cebeci achieve 

4 

accurate results at Reynolds numbers of order 10 in about 
0.9 s of CDC 6600 time. Since 1 s of CDC 7600 time equals 
about 5 s of CDC 6600 time and about 12 s of 1108 time, we 
see that the local finite-difference methods are comparable 
in speed to the global Chebyshev methods. However, it does 
seem that the local Chebyshev methods discussed above are 
faster by about an order of magnitude. But it should be 
pointed out that these comparisons are perhaps being done 
unfairly: M. Scott who wrote SUPORT informs us that SUPORT 
can be speeded up by about a factor 2 by simply changing 
the time-stepping algorithm; the Chebyshev codes can also 
be speeded up considerably in the areas of matrix set-up 
and full matrix pivoting. 



Table 5 . Some 

Results of the 

SUPORT Code 


Reynolds 

Error 

#of grid points 

#of iterations 

CDC7600 

number 




time/eigen- 

value 

o 

H 

1 

o 

H 

140 

4 

1.4 s 

10® 

1 

O 

H 

1400 

6 

23 s 
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1. NUMERICAL METHODS; EIGENVECTORS/ GROUP VELOCITY^ AND 
OPTIMIZATION METHODS 


Eigenvectors 

There are several efficient numerical methods to 
compute eigenvectors. In any of the standard finite-difference 
methods^ the eigenvector ip (y) is computed as part of the 
calculation of the eigenvalue X .-Also, in the local matrix 
iterative methods advocated in Sec. 6, the eigenvector is 
found as part of the iteration determining the eigenvalue X . 

A separate calculation of the eigenvector is required 
only for the case in which global matrix eigenvalue methods 
are used to determine X . In this case, one or two applications 
of the inverse iteration 


(A - XB)x 


(k+1) 


cBx 


(k) 


(7.1) 


where x^^^ is arbitrary usually determine the eigenvector 
to within roundoff error. The value of X used in (7.1) 
may be the calculated eigenvalue to machine precision — 
this normally does not cause any problems with singular 
matrices. For a 43 x 43 matrix, the calculation of the 
eigenvector by this inverse iteration scheme requires about 
75 ms using an assembly language fully pivoted LU algorithm 


on a CDC 7600. 


The computation of the eigenvector of the adjoint 
Orr- Sommer f eld equation may be done similarly (noting that 
the spectrum of the adjoint is the same as the spectrum of 
the Orr-Sommerfeld equation up to complex conjugation) . 


In practice, calculation of the adjoint egenf unction using 
the Chebyshev-spectral method requires about 150 ms using the 
assembly LU program because the matrix of the adjoint operator 
must be set up. 


Group Velocity 

The group velocity is of importance in relating the 
results of spatial and temporal stability theory eind in 
several optimization problems (see below) • In a layered 
flow with three-dimensional disturbances having wavevector 
( a , P ) and frequency o) ( a , g ) , the group velocity v^ is 


One way to compute the group velocity is pimply to compute 
the frequency o) for several nearby values of a , P and then 
use finite-difference approximations to v^. This procedure 
is neither efficient nor elegant. 

A much better way to compute the group velocity will 
now be described. It is usually faster than the crude finite- 
difference method described above (except possibly for the 
fastest local iterative eigenvalue solvers). We start by 
writing the Orr-Sommerfeld equation for three dimensional 
disturbances in the form 

= { (D^ - - g^) ^ - iR[ (au + gw - o)) (D^ - -. g^) 

« (au" + gw")]}i/^(y) = 0 (7.3) 

Taking the derivative of (7.3) with respect to a , we obtain 

= {4a(D^ - + iR(5 - |^) (D^ - 

- u" - 2a(au + gw -- uj)} ¥ (y) .. 
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Therefore, if x is the adjoint eigenfunction to tj» , we 
obtain 

3(1) ^ _ 4ig ^ /x{5(P^ - - 8^) - u" - 2a (gu + 3w - at) }\l> dy 

^ /X(D- r dy 

(7.5) 

There is a similar expression for 3w/33 • 

The computation of the group velocity using (7 ,5) 
requires little additional computational work to the 
calculation of the adjoint eigenfunction. 


Optimization Methods 

There are several kinds of information concerning 
the stability of a given flow that may be of interest. Some 
exaonples are; 

(i) A plot of a vs B at fixed R and frequency Re (a. 

(ii) Determination of that ot and 3 that maximize 
Im 0 ) at given R and Re o) . 

(iii) For given Im w a plot of a vs 3 at fixed R 

2 2 

or a + 3 vs R at fixed disturbance propagation 
angle or a neutral stability curve (a vs R for 
a two-dimensional disturbance with Im m = 0) . 

(iv) Determination of the critical Reynolds number, 

i.e. the smallest value of R at which there is a mode 
with Im 0 ) ==.0, 

(v) Computation of nonlinear ajid non-parallel flow terms. 

To illustrate how these problems can be solved efficiently 
on a computer, let us consider the solution of the problem (i) . 





An efficient procedure for obtaining the curve (i) 
of those a and B having Re w = f at a fixed Reynolds 
number R is to use Newton’s method as follows: 

(a) Starting with an approximation Cq, Bq to a point 
on the curve Re o) = f, we compute the group velocity 


at OL^, B^ and obtain the new approximation 

f - Re . 

g -' Rev 

|Re Vg! 

by shooting along the normal to rhe curve Re ci) = f . 


( 7 . 6 ) 


(b) Repeat step (a) until the approximation lies 
within a given error tolerance of the desired curve 
Re 0) = f. 

(c) Obtain an approximation to a new point on the 
curve Re 03 = f by shooting along the tangent to 

the curve from the previously found point; the tangent 
to the curve is in the direction Re(|^. " 

Repeat steps (a) and (b) • 


The solution of problem (ii) is obtained similarly. 

Here.it is possible to choose a new point on the curve 
Re 0 ) =5 f with a more favorable Im o) by use of either Newton’s 
method (which requires second derivatives of m ) or, say, 
cubic interpolation. One word of caution on this procedure 
for the solution of problem (ii)is that very small errors in 
the value of Re o) may confuse the search for the maximum of 
Im 03 • The reason is simply that in typical stability problems. 
Re 0) >> Im 0 ) . The moral is; the determination of optimal 

properties of the stability characteristics of flows requires 
extremely high accuracy calculations. 


81 


8. SUMMARY OF RESULTS 


A summary of our conclusions for a moderately fast 
stability code running on a CDC7600 computer is given 
in Table 6 for incompressible flow and Table 7 for 
compressible flow. 

Table 6. Typical computer timings for incompressible 

3 

Blasius flow at Reynolds number of order 10 with 
numerical errors of order lO”^ (Chebyshev polynomial 
methods with N = 34 polynomials) 

Operation Computer Time 


Temporal eigenvalue 


Matrix set-up 


0.09 

s 

Global search 


0.60 

s 

Local search 


0.04 

s 

Eigenvector 

Global search 
Local search 


0.03 

s 

Adjoint eigenvector 


0.09 

s 

Group velocity 


0.10 

s 

Nonlinear stability 

terms 

0.12 

s 

Nonparallel flow terms 

0.04 

s 

Spatial eigenvalue 

Global search 

Local search (first guess 

40 s 


by temporal 

code) 

0.4 s 

» 

Local search (guess known) 

0.04 

s 



Table 7. Typical computer timings for compressible flow 
stability calculations at Reynolds numbers of order 
10^ with errors of order 10~^ (Chebyshev polynomial 
methods with N = 34 polynomials) 


Operation 

Timing for 

Two-Dimensional 

Modes 

Expected Timing 
for Three- 
Dimensional Modes 

Temporal eigenvalue 

Global search 

30 s 

60 s 

Local search 

Direct 

2 s 

8 s 

Indirect 

0.25 s 

0.3 s 

Eigenvector and 

adjoint 

0.3 s 

0.35 s 

Nonlinear and 
non-parallel 
terms by transform 
method 

0.4 s 

0.5 s 


Final Remarks 

(A) Roundoff Error ; The effects of roundoff error 
become more acute as the matrix size increases. This 
effect is shown in Table 8. We conclude that high machine 
accuracy is necessary to even attempt calculations at large 
Reynolds numbers. 

(B) Accuracy of Profiles ; We have computed errors 
in the imaginary part of eigenvalues of as large as 10% 
with errors in the imposed profiles as small as .01%. 

However, these dangerous perturbations are of a very special 
kind wherein they are concentrated near the critical layer 

of the mode. (Obviously, if we perturb u(y) by an arbitrarily 
small amount, we may perturb u" (y)by an arbitrarily large 


Table 8 . Effect of roundoff error on the most unstable 

mode of plane Poiseuille flow for a = 1, R = 10000. 


M + 1 X (roundoff * lO"®) X (roundoff « 


20 

.23752685 

+ • 

•t. 00373451 

.23752676 

+ 

£.00373427 

23 

.23754139 

+ 

00383489 

.23752670 

+ 

£.00373982 

26 

.23749300 

+ 

£.00368897 

.23752646 

+ 

£.00373965 

38 

.23714159 

+ 

£.00352930 

.23752643 


£.00373966 

44 

.23348160 

+ 

£.00534311 

.23752643 

+ 

£.00373965 

50 

.23813295 

— 

£.00296263 

.23752655 

+ 

£.00373979 



amount at the same time. However, the dangerous perturbations 
cited above are perturbations in u" , but very noisy perturbations 
in u" near the critical layer.) 

In order to avoid these difficulties, it is necessary 
that the mean velocity profiles used in the stability 
calculations be very smooth and that the second derivatives 
of these profiles also be very smooth. 

(C) Accuracy of Kigenvalues ; Because the real parts 

of eigenvalues are typically much larger than the imaginai^ 
parts of the eigenvalues, it is necessary to maintain very 
high accuracy in stability calculations in order to get 
any meaningful result concerning instability. 

(D) Nonlinear and Non-Parallel Flow Terms ; The 
inclusion of these terms requires little additional work 
to that already expended in computing the eigenvalues of 
the flow. Therefore, it seems to make good sense to compute 
these terms . 

(E) Typical Computer Times ; It is realistic to 
expect a well-conceived stability code for incompressible 
flows to require less than 1/2 s per station of CDC 7600 
time; for three-dimensional modes of compressible flows, 
about 5 s per station is realistic. These estimates are 
for total computer time per station across a typical LFC 
wing provided either that an approximate eigenvalue is 
available at the leading edge or a previous station has 
been calculated nearby so that local methods may be used. 
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